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This proceeding is intended to be a first introduction to spectral methods. It is written around 
some simple problems that are solved explicitly and in details and that aim at demonstrating 
the power of those methods. The mathematical foundation of the spectral approximation is first 
introduced, based on the Gauss quadratures. The two usual basis of Legendre and Chebyshev 
polynomials are then presented. The next section is devoted to one dimensional equation solvers 
using only one domain. Three different methods are described. Techniques using several domains 
are shown in the last section of this paper and their various merits discussed. 



, I. INTRODUCTION 

o: . . . . . 

. This proceeding constitutes a short introduction to spectral methods. The aim is not to present an exhaustive 

■ mathematical presentation of those methods. Numerous books can be consulted for this purpose (see the bibliography 

O ^. for a sample of them) . The following material should instead be considered as a toolkit for implementing simple spectral 

D ' methods solvers. Thus, a particular emphasize will be put on practical examples. For the sake of simplicity, we will 

^0 ' restrict ourselves to one dimensional solvers. 

\^ Spectral methods are just one of the many ways to represent a function on a computer. The basic idea of all 
numerical techniques is to approximate any function by polynomials, those being the only functions than a computer 

. can exactly calculate. So a function u will be approximate hy u = X^^iLo "n^n where the $„ are polynomials and 

^ ' called the trial functions. Depending on the choice of trial functions, one can generate various classes of numerical 



o 

O 



techniques. For example, the finite difference schemes are obtained by choosing local polynomials of low degree. In 
all the following, we will be interested in spectral methods where the are global polynomials, typically Legendre of 
: Chebyshev ones. If spectral methods are basically more evolved than finite difference schemes, they have long prove 
' their ability to tackle a wide variety of problems. In particular, they allow to reach very good accuracy with only 
\^ , moderate computational resources. 

■ In the first part of this article, we present the basic properties of spectral methods. After a brief overview of 
' the foundations of the spectral expansion, we show the power of the method on simple but non-trivial examples. 
O^- In particular, we show that, for C°° functions, the error decays exponentially, as one increases the degree of the 
^ ' approximation. We then present some basic features of the two type of polynomials that are used in all the rest of 
(^JT), this article, Legendre and Chebyshev polynomials. 

^ The second section is devoted to the actual implementation of partial differential equation solvers. We restrict 

• ^ , ourselves to one dimensional equations on a bounded domain (i.e. x G [—1, 1]). Three different types of solvers are 
' presented: the Tau-method, the collocation method and the Galerkin method. They are explicitly implemented on a 
J-j _ simple example. 

. - - 1 The last part of this work is concerned with the implementation of multi-domain solvers. After having explained 
why this is often a valuable tool for the physicist, we present, once again, three different methods, all of them being 
tested an implemented on a simple, one dimensional equation, with a decomposition of space using two domains. 



II. FOUNDATIONS OF SPECTRAL METHODS 
A. Spectral expansion 

1. Orthogonal projection 

Let us consider an interval A — [a^min, a^max]- In order to talk about basis, one needs to define a scalar product on 
A. If w is a positive function on A, one can define the scalar product of two functions / and g, with respect to the 
measure w as being 

/ f{x)g{x)w{x)dx. (1) 

Using this scalar product, one can find a set of orthogonal polynomials p„, each of them of degree n. The set composed 
of those polynomials, up to a given degree A^ is a basis of Pat. 
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FIG. 1: u = cos^ {■Kx/2) — (a; + 1)* /8 and its projection of Chebyshev polynomials of degree smaller than 4 (left panel) and 8 
(right panel). 

One can then hope to represent any function u on A by its projection on the polynomials pn- Doing so, we define 
the projection of u simply by 

N 

PnU = ^ UnPn {x) , (2) 
n=0 

where the coefficients of the projections are given by u„ = ^^'^"-^ ^ xhe difference between u and its projection is 

{Pn,Pn) 

called the truncation error and one can show that it goes to zero when the order of the approximation increases : 

||u — Pnu\\ — > when A'' — > oo. (3) 

This convergence is illustrated on Fig. 1, where a test function u is plotted, with its projection for A'^ = 4 and N = 8. 
Chebyshev polynomials are used. Let us note that, even if u is not a polynomial function, no discrepancy can be seen 
by eye, for N as small as 8. 

This seems very appealing but for the fact that one needs to calculate the w„ by computing integrals like 

/ u (x) Pn (x) w (x) dx. If one needs to evaluate u at a lot of points in order to do so (like when using standard 
Ja 

Newton method), one would lose all the advantage of using spectral expansion. 



2. Gauss quadratures 

The solution to this problem is given by the Gauss quadratures. The theorem can be stated as follows. 
There exist N + 1 positive reals w„ and N + 1 reals A such that : 

f ^ 

V/ e P2N+S , f{x)w{x)dx = Y^f {Xn) Wn ■ (4) 

The Wn arc called the weights and the Xn the collocation points. The exact degree of applicability depends on the 

quadrature. The three usual choices are : 

• Gauss : S ~ I 

• Gauss- Radau : 5 = and xq = Xynin- 

• Gauss-Lobatto : S = —I and xq = .Tmin and xjv = Xy^n-x- 

Gauss quadrature is the best possible choice in terms of degree (it is not possible to find a quadrature with 6 > 1). 
However, the other quadratures have the property that the boundaries of the interval coincides with collocation points, 
which can be useful, to enforce boundary conditions for example. In all the following, we will use the Gauss-Lobatto 
quadrature. 
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FIG. 2: Same as Fig. 1 with also the interpolant of u. The collocation points are denoted by the circles. 

3. Interpolation 

If one applies the Gauss quadratures to approximate the coefficient of the expansion, one obtains 

N N 



= ;^X^w(a;j)Pn(a;j)wj with 7™ = (^j) "'j- (5) 

7n 



j=0 j=0 



Let us precise that this is not exact in the sense that u„ ^ u„. However, the computation of u only requires to 
evaluate u at the N + I collocation points. The interpolant of u is then defined as the following polynomial 



JV 



InU = UnPn {x) . (6) 



n=0 

The difference between /„u and P„u is called the aliasing error. The interpolant of u is the spectral approximate of 
u and one can show that it is the only polynomials of degree TV that coincides with u at each collocation point : 

[InU] {Xi) = U {Xi) Vi < N. (7) 

Figure 2 shows the same function as Fig.l but the interpolant is also plotted. One can see that, indeed Inu coincides 
with u at the collocation points that are indicated by the circles. Once again, even with as few points as A/^ = 8, no 
difference can be scxni bc;twc!c;n the various functions. 

Figure 3 shows the maximum difference between Jjvu and u on A, as a function of the degree of the approximation N. 
We can observe the very general feature of spectral methods that the error decreases exponentially, until one reaches 
the machine accuracy (here 10^^**, the computation being done in double precision). This very fast convergence 
explains why spectral methods are so efficient, especially compared to finite difference ones, where the error follows 
only a power-law in terms of N. We will later be more quantitative about the convergence properties of the spectral 
approximation. 

Let us note that a function u can be described either by its value at each collocation point u {xi) or by the coefficients 
of the interpolant Ui. If the values at collocation points are known one is working in the configuration space and in 
the coefficient space if u is given in terms of its coefficients. There is a bijection between the two descriptions and one 
simply goes from one space to another by using : 



1 ^ 

• Un= — Y2i'^{xj)pn (xj) Wj (configuration — > coefficient) 



N 

• u{Xn) = Y2^^Pi (coefficient — > configuration) 

Depending on the operation one has to perform, one choice of space is usually more suited than the other. For 
instance, let us assume that one wants to compute the derivative of u. This is easily done if u is known in the 
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FIG. 3: Maximum difference between Inu and w as a function of the degree of the approximation N. 



N=4 



N=8 





FIG. 4: First derivative of m = cos^ {Trx/2) — {x + 1)^ /8, interpolant of the derivative and derivative of the interpolant, for 
iV = 4 and AT = 8. 



coefRcient space. Indeed, one can simply approximate u' by the derivative of the interpolant : 



N 



U' (x) « [Inu]' = UnP'n i^) ■ 



(8) 



n=0 



Such an approximation only requires the knowledge of the coefficients of u and how the basis polynomials arc derived. 
Let us note that the obtained polynomial, even if it is a good approximation of u', is not the interpolant of u' . In 
other terms, the interpolation and the derivation are two operations that do not commute: (Inu)' =/= In (w'). This is 
clearly illustrated on Fig. 4 where the derivative of u = cos'^ {nx/2) — (x + 1) /8 is plotted along with the functions 
(Inu)' and In (m'). In particular, one can note that the functions used to represent u', i.e. (Inu)' does not coincide 
with u' at the collocation points. However, even with N = 8 only, the three functions can not be distinguished by 
eye. 

The maximum difference between u' and {Inu) , as a function of N, is shown on Fig. 5. Once, again, the convergence 
is exponential. 
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FIG. 5: Maximum difference between (Inu)' and u as a function of the degree of the approximation N. 

B. Usual families of polynomials 

In this section, we briefly present some basic properties of Legendre and Chebyshcv polynomials. Those two sets 
are the usual choice for non periodic problems. 

1. Legendre polynomials 

The Legendre polynomials, denoted by P„, constitute a family of orthogonal polynomials on [—1, 1] with a measure 
w = 1. The fact that the measure is so simple is one of the main advantage of working with Legendre polynomials, 
especially from the analytical point of view. 

The scalar product of two Pn is given by : 



The successive polynomials can be constructed by recurrence. Indeed given that Pq = 1 and Pi = x, all the P„ can 
be obtained by using 

(n + 1) P„+i (x) = (2n + 1) xPn (x) - nP„_i {x) . (10) 

It is then easy to see that the Legendre polynomials have the following simple properties : i) P„ has the same parity 
as n. ii) P„ is of degree n. iii) P„ (±1) = (±1)". iv) P„ has exactly n zero on [—1,1]. The first polynomials are 
plotted on Fig. 6. 

The values of the weights and collocation points can be written for the three usual quadratures : 

2 

• Legendre-Gauss : the nodes of Pm+i and Wi = o • 

il-xf)[P'^^,{x.f 

• Legendre-Gauss- Radau : xq = —1 and the Xi are the nodes of Pjv + Pjv+i- The weights are given by wq = 

2 , 1 

7T and Wi = TT. 

(N+lf {N+lf 

• Legendre-Gauss-Lobatto : a;o = — l,a;jv = — 1 and Xi are the nodes of PJ^. The weights are Wi = 

2 1 
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FIG. 6: The first Legendre polynomials, from Po to Pe- 



From the above formula, one can see that the position of the coUocation points are not analytical and they have to 
be computed numerically, which is one of the main shortcoming of Legendre polynomials. 

One can also derive the action of some linear operation in the coefficient space. Consider a function / given by its 
coefficients : / = Y^n=Q '^nPn (x) and _ff be a linear operator acting on / so that Hf = X^^_o bnPn {x)- The relation 
between the a„ and fe„ can be explicitly written in some cases. For example : 

• If -ff is the multiplication by x then : 

• If if is the derivation : 

N 

6„ = (2n + l) '^p- (12) 

p=n+l,p+n odd 

• If is the second derivation : 

N 

6„ = (n+l/2) \p{p+l)-n{n+l)]ap. (13) 

p— n4-2,p4-neven 

Those kind of relations enable to represent the action of if as a matrix acting on the vector of the a„, the product 
being the coefficients of if/, i.e. the 6„. 



2. Chebyshev polynomials 



The Chebyshev polynomials T„ are an orthogonal set on [—1, 1] for the measure w = ^ More precisely one 

V 1 — 

has 



' ^^dx=^{l + 6on)Smn- (M) 

-1 vl — a;^ ^ 
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FIG. 7: The first Chebyshev polynomials, from To to Te. 



Chebyshev polynomials can be computed by knowing that Tq = 1, Ti = x and by making use of the recurrence : 

T„+i {x) = 2xTn (x) - T„_i (x) . (15) 

It follows that i) Tn has the same parity as n. ii) T„ is of degree n. iii) T„ (±1) = (±1)". iv) T„ has exactly n zero 
on [—1. 1]. The first polynomials arc plotted on Fig. 7. 

The weights and collocation points associated with Chebyshev polynomials can be computed : 

• Chebyshev-Gauss : Xi = cos and Wi = ^ 



2N + 2 ' N + 1 

Chebyshev-Gauss- Radau : Xi = cos . The weights are wq = -r-r^ — - and 



2Ar-hl ° 2N+1 2N+1 

• Chebyshev-Gauss-Lobatto : Xj = cos—. The weights are wq = wn = and Wi = —. 

Contrary to Legendre polynomials, the position of the collocation points are completely analytical which can somewhat 
simplify the computational task. In all the following, Chebyshev polynomials are used, except when otherwise stated. 

Once again, one can find the relation between the coefficients o„ of a function / and the coefficients 6„ of Hf, 
where is a linear operator. For example : 



If H is the multiplication by x then 



If H is the derivation : 



bn = ^[{'i- + Son-i)an-i+an+i] (n > 1) . (16) 



2 ^ 

[^ + Oon) , , , 

p=n-\-l,p-\-nodd 



If H is the second derivation : 



N 



'■"(TT^) Jl ^(v^-n'W (18) 

^ ^ p— n+2,p+n even 
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C. Convergence properties 



As already seen, for sufficiently regular functions, the difference between u and its interpolant I^u goes to zero 
exponentially. This statement can be made more precise. Let us consider a C™ function u. Upper bounds on the 
difference between u and I^u can be found for various norms and choice of polynomials. 

• For Legendre : 



\\Inu — u\\j^2 < 



C 



jym-l/2 



El 

fe=0 



,(fe) 



For Chebyshev : 



||7jvw- mI 



< 



C 



— T 



fe=0 



hi. 



(19) 



(20) 



< 



C 



-1/2 ^ I 
fe=0 



L2 



(21) 



Without going into to much details, one can note that the errors decay faster than some power of A'^. The upper 
bounds decay faster for more regular function (i.e. for higher m). An interesting case is the one of a C°° function. In 
such a case, the error decays faster than any power of A'' and thus like an exponential. One talks of an evanescent 
error and this is the case previously observed. 

A limit case can be provided by a discontinuous function. Consider a step function u such that u (a; < 0) = 1 and 
M (a; > 0) = 0. If one tries to interpolate this function the convergence theorems previously exposed can not ensure 
convergence. This reflects on Fig. 8 where the step function is presented with interpolants for various number of 
points. As the number of collocation points increases, the maximum difference between u and Ijqu stays constant (the 
amplitude of the oscillations can be seen to be constant). However, more and more oscillations are present, as the 
Chebyshev polynomials do their best to approximate the discontinuous step function. In a sense, the approximation 
is still better and better, the support of the difference being smaller and smaller. In other words, the integrated error 



I dx still goes to zero when N increases, as can be seen on Fig. 9. However, this convergence follows 



only a very slow-decaying power-law. 

The simple example of the step hmction is a very general feature of the so-called Gibbs phenomenon, which occurs 
anytime one tries to interpolate a function that is not C°°. When this happens, the error is no longer evanescent and 
only converges as a power-law. In order to maintain evanescent errors, one should try to avoid any Gibbs phenomenon. 
In some cases, as will be seen later, this can be achieved by using a multi-domain decomposition of space. 



III. DIFFERENTIAL EQUATION SOLVERS 
A. The weighted residual method 



The type of problems with are concerned with in the section are ordinary differential equations, in one dimension 
only, on a bounded domain at the boundaries of which some conditions are enforced on the solution. Mathematically, 
one considers the following system : 

Lu (x) ^ S (x) for xeU (22) 
Bu {y) = ioT yedU (23) 

where L and B are linear differential operators. A function u is then an admissible solution of this system, if and only 
if i) it satisfies Eq. (23) "exactly" (i.e. to machine accuracy) ii) it makes the residual R = Lu — S small. In order to 
quantify what this "small" means, the weighted residual method relies on A'' + 1 tests functions ^„ and one asks that 
the scalar product of R with those functions is exactly zero : 

(a,i?) = 0, \/k<N. (24) 

Of course as N increases the obtained solution is closer and closer to the real one. Depending on the choice of spectral 
basis and of test functions, one can generate various different types of spectral solvers. In the following, three of the 
most used ones are presented and applied to a simple case. 
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FIG. 8: Step function and some iutcrpolants for N = 7, N = 15 and N = 31. As N increases the maximum difference stays 
constant and the number of oscillations increases. 
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FIG. 9: J^^ \Inu — m| da; as a function of N, u being a step function. The convergence obeys a power-law which decays slower 
than N~^. 



We propose to solve the equation : 



B. A test problem 



-r^ - 4— + 4^ = exp (x) + C, 

ax^ ax 



(25) 



with X e [—1, 1] and C = — 4e/ (l + e^). As boundary conditions, we simply ask that the solution is zero at the 
boundaries : 



m(-1)=0 and w(l) = 0. 



(26) 
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Under those conditions, the solution is unique and analytical : 

, , , , sinh (1) C 

Wsoi {x) = exp (x) - exp (2a;) + -. 



(27) 



Let us note that this solution is not a polynomial. 

d^ d 

With our notations the linear operator L is -p-^ — 4- — |-4Id. Using the elementary linear operations seen in Sec. II B, 

da;^ da; 

one can construct the matrix representation of L, which will come handy in the implementation of the various solvers. 

N N N 

Let us recall that if u = UiTi (x) then Lu — LijUjTi (x). 

i=0 i=0 j=o 

In this particular case, and for N = 4, one finds : 



/4 -4 4 -12 32 \ 
4 -16 24 -32 
4 -24 48 
4 -32 

\0 4 / 



(28) 



The Tau-method 



In this method, the test functions ^„ are chosen to be the same as the spectral functions of decomposition. Let us 
assume that one uses Chebyshev polynomials Tj. The residual equations (24) are then : 

{T„,Lu-S) = \/n<N. (29) 

If one uses the definition of the matrix Lij , this set of equations can be written as 



N 

E 

3=0 



Lr. 



iUj Sr. 



yn<N 



(30) 



where the s„ are the spectral coefficients of the source S. 

However, due to the presence of homogeneous solutions of L, this set of A'^ + 1 equations is degenerate and one 
must impose the boundary conditions before solving it. In the Tau-method, the boundary conditions are enforced as 
additional equations. In our particular case they can be written in a very straightforward manner : 



u{x = -l) = 



u(a; = -1-1) = 



N 

i=o 

N 
.5=0 



(31) 



One has then to relax the two last residual equations and to replace them by the two boundary conditions in order 

to get an invcrtiblc system for which the unknowns arc the w,„. Relaxing the last two equations is not an issue. 
Indeed, if the functions are regular enough, the coefficients are rapidly decreasing and thus the obtained solution 
should converge nicely to the exact solution. 

In our example, and for A'' = 4 the system reads as : 



(32) 



Once inverted, the coefficients of the solution can be found mq ~ 0.146; ui ~ 0.079; W2 — — 0.122; M3 ~ — 0.079; U4 ~ 

-0.024 (for N = 4). 

The numerical solution obtained is compared to the exact one on Fig. 10. If the difference is easily seen for = 4, 
even with as few points as A'' = 8, there is not difference noticeable by eye between the two curves. One can also 
note, as expected when using the weighted residual method, that, for all N, the numerical solution "exactly" fulfills 
the boundary conditions. 
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FIG. 10: Exact solution (solid line) and numerical one (dashed line and cross symbols), for the Tau-method, for AT = 4 and 
N = 8. 



D. The collocation method 



In this method, the test functions are no longer the functions of the spectral expansion, but are chosen to be zero 
but at one collocation point = ^ {x — Xn)- The residual equations (24) are then equivalent to : 

Lu (xn) = S (xn) Vn < N. (33) 

Given the expression of the coefficients of Lu in terms of the matrix i^, one can write those equations as a system 
for the unknowns Un '■ 

N N 

J2 ^ioTi {Xn) Uj = S {Xn) ^Tl < N. (34) 

As for the Tau-method, this system is degenerate and must be supplied with the boundary conditions. One again, 
they are enforced by two additional equations, which are the same than for the Tau-method (i.e. Eqs (31)). However, 
instead of relaxing the last two equations of (34), one has to relax the first and the last ones. Indeed, those are the 
equations that involved the extremal points, at which the boundary conditions have to be fulfilled. 
For N = 4, the system for the collocation method is : 



/I 
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-1 


1 \ 
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(35) 



and the coefficients of the solution are : uq ~ 0.188; wi ~ 0.089; M2 — — 0.157; M3 ~ — 0.089; M4 ~ —0.031. 

Once again, the numerical solution and the exact one are shown on Fig. 11 and the same properties than for the 
Tau-method can be observed. 



E. The Galerkin method 



The basic idea of the Galerkin method is to expand the solution, not in terms of usual orthogonal polynomials, but 
on some linear combinations of polynomials that fulfill the boundary conditions. One then talks of Galerkin basis. 
The particular choice of basis is of course important and it is rather hard to give a general recipe. However, it is 
usually better if the Galerkin basis can be easily written in terms of the original basis. 

For the boundary conditions of our example (26), it is easy to see that the following choice of Galerkin basis Gj is 
a valid one : 

• G2k (X) = T2k+2 (x) - To (x) 
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FIG. 11: Exact solution (solid line) and numerical one (dashed line and cross symbols), for the collocation method, for iV = 4 
and N = 8. 



• G2k+i {x) = T2fe+3 (x) - Ti {x). 



Let us note that, in order to maintain a consistent degree of approximation, one has to consider only N —1 Galerkin 
polynomials (so that the Chcbyshcv polynomial of higher degree remains Tj\f). This is evident on our example but is 
a general feature. The reader should also note that the Gi, in general, are not orthogonal polynomials. 

It is convenient to define the transformation matrix as the matrix that relates the Gi to the standard polynomials, 
let us say the Tj. This transformation matrix M is not square but a (A'' + 1) x (A'' — 1) matrix and is defined by 



N 
3=0 



yi<N -2. 



(36) 



In our example M is easily constructed and for iV = 4 one gets 



Mi, 



/-I -1\ 
0-10 

1 
1 

Vo 1/ 



(37) 



N-2 



As already mentioned, the solution u is seek in terms of the Galerkin basis : u = u^Gk (x). The transformation 

fe=o 

matrix enables us to get u in terms of Ti, and then to use the matrix Lij to get the expression of 

N-2 N N 



Lu ix)^J2''kJ2J2 MjkLijTi (x) . 

=0 j=0 



(38) 



fe=0 



In order to use the weighted residual method, one has to precise our choice of test functions. In a standard Galerkin 
method, the ^„ are chosen to be the Galerkin coefficients so that the residual equation reads ((?„, R) < N — 2. 

{Lu, Gn) can be computed by using Eq. (38) and by expressing G„ in terms of Chebyshev polynomials, using, once 
again, the transformation matrix. 

The situation for the source term is a bit different. Indeed, the source S has no reason at all to verify the same 

JV 

boundary conditions as the solution u. So, one has to expand S on the standard basis : S = ^^SiTj. The scalar 

i=0 

product {S, Gn) can then be obtained by using M to go from the Gi basis to the Tj one. 
Putting all the pieces together, the Galerkin system then reads : 



N-2 N N N 

"fe MinMjkLij {Ti\Ti) = ^ MinSi {Ti\Ti) , Vn < TV - 2. 

fe=0 i=0 j=0 i=0 



(39) 
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FIG. 12: Exact solution (solid line) and numerical one (dashed line and cross symbols), for the Galerkin method, for N = 4 
and N = 8. 



This system is well-posed and can be inverted to get the unknowns, which are the coefficient of u on the Galerkin 
basis, i.e. the uf. Once those are known, the transformation matrix can be used, one last time, to get the solution 
in terms of the standard spectral basis : 



N /N-2 



"(^) = E E^»^- M- (40) 

i=0 \n=0 / 

In the case of our test problem, for A'' = 4 the system is : 

'27r -47r -Air 

Stt -Stt I I I = I -1-705 I . (41) 
Stt -26n^ 

The Galerkin coefficients arc Uq' ~ — 0.160; ~ — 0.092; ~ —0.029 and the standard ones ~ 0.189; Ui ~ 
0.092; U2 — —0.160; U3 ~ —0.092; U4 ~ —0.029. The numerical solution is compared to the exact one on Fig. 12. 
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F. The methods are optimal 



On Fig. 13, the maximum error between the numerical solution and the analytical one is shown, as a function of 
the number of coefficients N . The three methods exhibit similar features : an exponential decay before reaching the 
machine accuracy and being limited by round-off errors. All the functions involved in our test being C°°, it is not 
surprising that the error is evanescent. Also plotted on Fig. 13 is the maximum difference between the exact solution 
and its interpolant. The behavior is very similar to the three other curves. Thus, it seems that all the three resolution 
methods do not introduce more error than the one that would be already done when interpolating the exact solution. 

This notion can be made more precise and is called optimality. Let us assume that the exact solution of an equation, 
Wexact, is known. One can compared it to its interpolant /jv^exact and to the obtained numerical solution Wnum- The 
resolution method is then said to be optimal if and only if max(|ucxact — ^wUcxactl) and max(|?ioxact — u,ium.|) have 
the same behavior when N ^ 00. Of course, this is usually a difficult quality to assess, given that the exact solution 
is rarely known. However, on our particular simple example, all the three methods presented are optimal. 



IV. MULTI-DOMAIN SOLVERS 



A. Motivation 



As previously seen, spectral methods loose some of their appeal when they have to deal with non C°° functions. 
Indeed, due to the Gibbs phenomenon, the convergence of the numerical solution to the real one is drastically slowed 
down. 
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FIG. 13: Maximum relative difFereuee between the exact solution and the numerical one as a function of the number of 
coefficients, for the Tau-method (dashed line), the collocation one (dotted one) and the Galerkin method (dashed-dotted line). 
Is also shown, the relative difference between the exact solution and its interpolant. 



Nevertheless, they are a great deal of physical problems where such discontinuities occur. We can for example 
mention the shock front in fluid dynamics or the surface of a strange star where a jump of matter density occurs. 
Prom a mathematical point of view, it is also sometime useful to use different variables and functions in various regions 
of space, thus creating numerical discontinuities. 

For all those reason, one has to introduce several numerical domains to cover the physical space of interest. If one 
can ensure that the various discontinuities lie at the boundaries, then, in each domain, one will deal only with C°° 
functions, thus recovering exponential convergence. 



B. Test problem 



As in Sec. Ill, the various multi-domain solvers are presented on a simple test problem. The physical space is the 
interval — 1 < x < 1. On this interval, one considers the following equation : 

-g^+4« = 5 (42) 

with the boundary conditions u{—l) = and u{l) = 0. The source 5 is a step function : 5 (a; < 0) = 1 and 
S (x > 0) = 0. Under those conditions, the analytical solution is given by : 

u{x<Q) = i - + S-e^^ exp (2a;) + exp (-2a;) (43) 

M (x > 0) = B+ ^exp (-2a;) - ^ exp (2a;)^ , (44) 

where the constants are 

1 e2 



B- 



;(l + e2) 8(l + e4) 



/ 1 
B+ = — ' 



,(l + e4) (l + e2). 
As for the one-domain case, this solution is not polynomial. 
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Number of coefficients 



FIG. 14: Error done when using a Tau-nicthod with a discoutinuous source, as a function of the number of coefficients. The 
convergence is no longer exponential but follows a power-law in . 

Of course, due to the discontinuity of the source, the solution is only at x — 0. This causes a great lose of 
accuracy for the one-domain solvers presented in Sec. III. For example, the error made when solving Eq. (42) with 
the standard Tau-method is presented on Fig. 14. As expected, even if one converges to the right solution, the decay 
of the error is no longer exponential but only follows a power-law (like N~'^ in this particular case). For a practical 
point of view, the lose of accuracy is spectacular : with as many coefficients as = 100 one does not oven reach a 
precision of 10"'^. This is to be compared with the results of Sec. Ill, where machine precision could be reached with 
only 25 coefficients (see Fig. 13). 

In order to cope with this problem, one can introduce two subdomains of the interval : 

• Domain 1 : a; < described by the variable xi =2x + 1 (— 1 < xi < 1). 

• Domain 2 : a; > described by 2:2 = 2a; — 1 (—1 < X2 < 1). 

In each subdomain, the functions arc expanded on the basis polynomials, with respect to the auxiliary variable. For 
instance, in domain 1, the solution is given by : 

N 

u{x<0) = Y, ^iTi (xi (x)) . (45) 

When implementing various operators, one has to remember that one is working with the auxiliary variables. For 

, d „ d 
example — = 2— — . 
dx dxi 

C. A multi-domain Tau-method 

In each domain, one writes the usual residual equations for the Tau-method (see III C). For instance, in the domain 
1, this gives raise to 

N 

(T„, R)=Q=^Y, Lniu] = si (46) 

i=0 

This is a set of A'' -|- 1 equations among which one has to relax the last two, in order to impose various boundary and 
continuity conditions. L is the differential operator expressed with respect to Xi. The same is done in the second 
domain. 
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As for the usual Tau method, the boundary conditions are enforced by two additional equations. In order to close 
the problem, one has to impose that the solution and its first derivative with respect to x are continuous at the 
boundary between the two domains : 

• {xi = l,x = 0) = {x2 = — 1, a; = 0) 

._(^, = 1,^ = 0) = — (X2 = -I,x = 0) 

We are left with a system of : 2A'' — 2 residual equations, 2 boundary conditions and 2 matching conditions, for a 
total of 2 A'' + 2 equations. The unknowns are the coefScients of the solution in both domains which are in the right 
number 2N + 2. For our test problem, when using A' = 4 in both domains, the system to be inverted is : 
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(47) 



The solution is then given in both domains by 

• ~ 0.082; ~ 0.044; ~ -0.036; {t^ ~ 0.0018; ul ~ -0.0007 

• ~ 0.038; {tf ~ -0.044; ui ~ 0.008; «i ~ -0.0018; u| ~ 0.00017 
The convergence properties of this method will be presented in Sec. IV F. 



D. Homogeneous solution method 



This method is the closest to the analytical way of solving ordinary differential equations in various domains. The 
idea is just to compute some particular solution in both domains, to get the homogeneous solutions and to do the 
appropriate linear combination to fulfill the boundary conditions and continuity requirements. 

Most of the time, the homogeneous solution are known analytically. Gcnerically, their number is given by the 
order of the operator L of the equation. In our particular example, L is of second order and one gets 2 homogeneous 
solutions in each domain. One has to compute the coefficients of those solutions with respect to the auxiliary variables 
x\ and X2- In our case, one finds : 

• hi {x) = exp (2a;) = exp (a;i — 1) = exp {x2 + 1) 

• /i2 {x) = exp (—2a;) = exp {—xi + 1) = exp {—X2 — 1). 

For instance, for hi this gives raise to the following coefficients, in both domains : 

• in the domain 1 : (0.47; 0.42; 0.1; 0.017; 0.0020) 

• in the domain 2 : (3.4; 3.1; 0.74; 0.12; 0.015) 

In order to get one particular solution in each domain, one can use whatever method already presented in the one- 
domain case (sec Sec. III). The conditions enforced at the boundaries of each domain have no importance because 
they will be modified by linear combination with the homogeneous solutions. For instance, using a Tau- method in 
each domain, and demanding that the first two coefficients of the particular solutions are zero (equivalent to imposing 
two boundary conditions) gives the following system in the first domain (which reduces to a (A^ — 1) x (A^ — 1) system) 



/ 1 \ 
1 ' 
4 -16 -128 

4 -96 
VO 4 -192/ 



/~9o\ 
91 

92 

93 

\~9aJ 




1 


Vo/ 



(48) 
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which enables us to obtain the following particular solution, in domain 1 = 0; = 0; 52 — —0.53; ^3 = 0; 54 = 
-0.001. 

So, in each domain, the solution is given by the sum of the particular solution g and the two homogeneous ones hi 
and /i2, with unknowns coefficients a* and /3% 

u (x) = {xi) + a'h\ {Xi) + {xi) (49) 

where i denotes the domain. 

In order to determine the appropriate vahies of the 4 coefficients a' and one has to impose : 

• Boundary condition at a; = — 1 i.e. at xi = — 1 

g^ {-l) + a^hl {-l) + p^hl{-l)=0 

• Boundary condition at a; = 1 i.e. at 0:2 = 1 

• Matching of the solution at x = 0, i.e. at xi = 1 and X2 = — 1 

g^ (1) + a'hl (1) + P^hl (1) = g^ (-1) + a^hj (-1) + P^hj (-1) 

• Matching of the first derivative at a; = 0, i.e. at xi = 1 and a;2 = — 1 

g'^ {l)+a^K^ {1)+I3^h'^' (1) = 5" {-l) + a^hf (-1) + P^h'i (-1) 



If the problem is well-posed, this system of 4 equations admits a unique solution for the 4 coefficients of the 
homogeneous solutions, thus allowing us to construct the solution of the equation in all space. For instance, in our 
particular case, for N = A, the matching system is 



/7.39 0.135 \ 
1 1 -1 -1 1 I 

1 -1 -1 1/32 
\ 0.135 7.39/ 

which admits the solution : 0^ = 0.0048; = 0.14; /J^ = 0.094; 

solution, in each domains are then : 




(50) 



= —0.0017. The coefficients of the numerical 



ul ~ 0.083; u\ ~ 0.044; ~ -0.036; ul ~ 0.0018; 



-0.0008 



• r4 - 0.038; uf ~ -0.044; ~ 0.008; ui ~ -0.0018; uf ~ 0.00016 
Numerical accuracy of this method will be shown in Sec. IV F. 



E. Variational formulation 



Contrary to all the methods previously presented, the variational one is only easily applicable when one is using 
Legendre polynomials. Indeed, as will be seen later, it requires that the measure is w (x) = 1. 
The starting point is the usual residual equation, written here explicitly in its integral form 



(^, i?) = =^ j ^ {-u" + 4u) da; = y" ^Sdx. 
The term involving u" is then integrated by parts, which gives : 

[-4"'] + / ^'u'dx + 4 / ^udx = / ^Sdx. 



(51) 



(52) 



The method is easier to implement when, as for the collocation method, the test functions are : ^„ = i5 (a; — a;„), 
being zero but at one collocation point. 
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In the following, it will also be eonvenient to define the derivative operator D in the configuration space. It relates 
the values of the function at each collocation point to the ones of its first derivative via : 



N 



9' {xk) = ^Dkjg{xj). 
3=0 

The matrix associated to D is easily constructed. For example, for A'' = 4 and Legendre polynomials, one gets : 

/ -10 13.5 -5.33 2.82 -1 

-2.48 3.49 -1.53 0.52 
0.75 -2.67 2.67 -0.75 
-0.52 1.53 -3.49 2.48 
V 1 -2.82 5.33 -13.5 10 / 



(53) 



(54) 



Using this operator and the Gauss quadrature rule (4), one can compute the integrals involved in Eq. (52) : 

N N N 



j i'nU'dx = ^ {Xi) u' {Xi) Wi = ^^ DijDinWiU {Xj) 



i=0 3=0 

inW^X = ^ {Xi) u{Xi)Wi = U {Xn) Wn 



i=0 

N 



i=0 

N 



/ ^nSdx = {xi) S {Xi) Wi = S (Xn) Wn- 



(55) 
(56) 
(57) 



For the points strictly inside each domain, the integrated term [— ^m'] = so that the residual equations take the 
following form : 



N N 



DijDinWiU (Xj) + 4u {Xn) W„ = S {Xn) Wn, < n< N 

i=o i=o 



(58) 



which is a set of 2N — 2 equations. 

The situation of the point lying at the boundary between the two domains is a bit tricky. Indeed the same physical 
point can be described either : 

• as the last point of domain 1 : n = N and [— ^m'] = —u'^ {xi = 1; a; = 0) so that the residual equation is : 

TV TV 

u'^ (a;i = 1) = ^ ^ DijDiNWiv} {xj) + 4m^ (a;jv) wn - {xn) wn- (59) 

i=0 j=0 

• as the first point of domain 2 : n = and = u'^ {x2 = — 1; a; = 0) and the residual equation is then 

N N 

u'^ (a;2 = -1) = - ^ ^ DijDioWiV? (xj) - 4u^ (xq) wq + (xq) wq. (60) 

i=0 j=0 

This duality can be used to impose implicitly the continuity of the first derivative of u : 
v!^ (xi = 1; a; = 0) = {xi = -1; a; = 0) =4- 

N N N N 

DijDiNWiU^ (xj) + 4u^ (a; at) wat + ^ ^ DijDioWiU^ (xj) + Au^ (xq) wq 

i=0 j=0 i=0 3=0 

= S'^ {xn)wn + {xo)wo. (61) 

The other equations are i) the boundary condition on the left {xo) = ii) the boundary condition on the right 
V? {xn) =0 and iii) the matching of the solution at the boundary between the two domains v} {xn) = u"^ {xq)- This 
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is a total number of 2N + 2 equations for the 2N + 2 unknowns, which in this case, are the value of the solution at 
each collocation points, the (a;„). 

For our example and A'' = 4 in each domain, the system is : 
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Once it is inverted, the solution is found in the configuration space : 

• In domain 1 : it^ (.t„) = (0; 0.06; 0.12; 0.12; 0.092) 

• In domain 2 : (a:„) = (0.092; 0.064; 0.030; 0.0090; 0). 

As for the other methods, the convergence properties are shown in the next section. 

Before finishing with the variational method, it may be worthwhile to explain why Legendre polynomials were used. 

Suppose one wants to use Chebyshev polynomials. The measure is then w {x) = . When one integrates the 

vl — 

term in u" by part one then gets 

j -u"fwdx = [-u'fw] + j u'fw'dx (63) 

Because the measure is divergent at the boundaries, it is difficult, if not impossible, to isolate the term in u'. On the 
other hand this is precisely the term that is needed to impose the appropriate matching of the solution. There might 
be ways around this but this explains why the variational method has been presented with Legendre polynomials. 

F. Discussion 

On Fig. 15, we show the maximum relative difference between the numerical solution and the analytical one. 
Contrary to the case of just one domain (see Fig. 14), exponential decay of the error is recovered, as expected, and 
machine accuracy is very rapidly reached. This illustrates the usefulness of a multi-domain decomposition. 

From a numerical point of view, the method based on the homogeneous solutions is different from the two others. 
Indeed, when using it, the problem is treated domain by domain. The number of system to be solved is equal to the 
number of domains but each of those systems is roughly of size N"^ only. This is to be compared to the Tau-method 
and the variational one where only one global system is inverted but which is of size {NNr,y. So, if the number of 
domain is important, one may consider the homogeneous solution method instead of the two others, which would give 
rise to huge systems. 

On the other hand, the Tau-method and the variational one, do not require the knowledge of the homogeneous 
solutions which is a good point, especially when the differential operator L is complicated. The variational method 
may seem a little more involved, especially as it requires to work with Legendre polynomials but, from a mathematical 
point of view, it is the only one that has been prove to be optimal. 

V. CONCLUSION 

Our introductory tour of spectral methods is now over. After having presented the mathematical foundations of 
this class of method, along with two families of standard polynomials, Legendre and Chebyshev ones, two simple test 
problems were treated by means of various methods. Three ways of solving an ordinary differential equation on a 
single domain were presented and tested : the Tau method, the collocation one and the Galerkin method. In the 
last part, a discontinuous source was treated by means of a multi-domain decomposition. Again, three methods were 
implemented and discussed : a multi-domain tau method, one based on the knowledge of homogeneous solutions and 
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FIG. 15: Relative differ ence between the analytical solution and the numerical one, for the three multi-domain methods exposed. 
The error is evanescent. 

a variational method using Lcgendre polynomials. This simple examples were intended to illustrate the usefulness of 
spectral methods to reach very good accuracy with moderate computational resources. 

However, this work is far from covering the huge field of all spectral methods. Among the many aspects that were 
ignored one can mention the use of Fourier transform for periodic functions. We also have left out all issues concerning 
compactification of space, by means of auxiliary variables like 1/r. Three dimensional problems were also left out as 
were non-linear problems. Finally, we did not present any problem involving time evolution. Should the reader feel 
the urge to learn more about all this, we recommend that he consults some of the books given in the bibliography. 
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